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ABSTRACT 

Approximate analytical formula for density distribution in differentially rotating stars is de- 
rived. Any barotropic EOS and conservative rotation law can be handled with use of this 
method for wide range of differential rotation strength. Results are in good qualitative agree- 
ment with comparison to the other methods. Some applications are suggested and possible 
improvements of the formula are discussed. 
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1 INTRODUCTION 

Theory of self-gravitating rotating bodies seems to be an unlimited 
reservoir of difficult problems hardly tractable even under severe 
simplifications. It is a subject of scientific effort since 1742 when 
Maclaurin has initiated this field by his studies on incompressible 
rotating ellipsoids I Maclaurinl ll742h . Development of modern nu- 
merical calculations resulted in progress in practical applications 
nowadays such as e.g. 3D hydrodynamical simulations of rotation 
of complex objects. 

Analytical approach has succeeded for constant density, 
incompressible bodies. Work of Maclaurin, Jacobi, Poincare, 
Schwarzschild and many others has explai ned th e behaviour of 
those objects almost completely lLvttletor3 ll953l). Behaviour of 
slowly rotating polytropes has been calculated by Chandrasekhar 
1 1936). By applying the differential equation of hydrostatical equi- 
librium modified by rotation he reduced the problem to an ordi- 
nary differential equation. This method however works only for a 
uniform rotation. This list would be incomplete without the Roche 
model. It's simplicity makes it a very powerful tool for understand- 
ing behaviour of rotating objects. Present computational methods 
allow one to handled numerically two- and three-dimensional prob- 
lems with complicated governing equations 

In this paper we present simple analytical approach which can 
treat differentially rotating compressible barotropic stars in case of 
slow or moderately fast rotation. This model could fill a gap be- 
tween simple analytical methods used for e.g. Maclaurin spheroids 
or Roche model, and complicated numerical methods such as e.g. 
HSCF lHachisull986[). or those applying str aightforward Newton- 
Raphson technique lEriguchi & Mullerll985l) . 



2 FORMULATION OF THE PROBLEM 

We attempt to find a density distribution (iso-density contours) of 
a single self-gravitating object under the following assumptions: 

(i) Barotropic EOS p = p(p) 

(ii) Simple rotation v = r f2 e^, with angular velocity dependent 
only on the distance from rotation axis — fi(r) 

(iii) Newtonian gravity 

(iv) Axisymmetric density distribution 

(v) We seek solutions for stationary objects in full mechanical 
equilibrium, i.e. all quantities are time-independent 

With properties (i)-(v) satisfied, the Euler equation becomes, 
in cylindrical coordinates (r, z, (f>): 



r n(r, zf e r = - Vp + V$„ 
P 



(1) 



Continuity equation is then fulfilled automatically. Introducing cen- 
trifugal potential: 



* c (r) 



and enthalpy: 



Hp) 



f Q(r) 2 df 



(Up 

J p 



we get a simple equation: 

V[ft(p)+* s + $c] 

with a solution 

C-- 



const 



(2) 



(3) 



(4) 



(5) 
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Hp) + * s + * 

Equation J5J is the most important equation in the study of the 
structure of rotating stars under conditions (i)-(v). We define the 
integration constant in eq. {3J to be such that the enthalpy satisfies 
the condition h(p — 0) = 0. The only term which we haven't 
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specified yet is the gravitational potential $ 9 . If we use the Poisson 
equation: 

A$ 9 =4 7rGp (6) 

where G is the gravitational constant, the equation J5J becomes 
a non-linear second-order differential equation. This form, how- 
ever, is very inconvenient, because we have to specify boundary 
conditions at a surface of the star, 1 which is unknown a priori. 
More powerful is an integral form of eq. J5J obtained by substitu- 
tion: 

This integral form has been used in very successful num erical 
algorithm developed originally by Ostriker & Mark 1 1968), and 
recently improved by Hachisu ( 1986) and by Eriguchi & Miiller 
( 1985). This form will be also used to derive our approximation 
formula in the next section. 



3 APPROXIMATION FOR THE DENSITY 
DISTRIBUTION 

The integral equation form of eq. J5J is: 

h{p) + TZ(p) + $ c = C (8) 

where 1Z is the integral operator acting on the density performing 
the integration on the right-hand side of eq. J7) over entire volume 
of the star. We define the surface of the star to be manifold consist- 
ing of points where p = 0, h = 0. Explicit form of the operator TL 
in terms of coordinates will not be needed. For a given EOS and for 
a fixed rotation law i.e. for given functions h(p) and <£> c , the only 
free parameter is the constant C. The values of C label a family of 
the stellar models with the same EOS and rotation pattern, which 
differ in total mass and maximum density 2 etc. 

Eq. (8j has a form of t he Hammerstein non-linear integral 
equation iHammersteidl 1 93cT) and can be rewritten in a canonical 
form: 

f = K [F(f)] (9) 

where: 

f = C-$ c -h(p), F(f) = h- 1 (f + $ c -C) (10) 

In case of linear function F, eq. {5) could be easily solved by the 
von Neumann series. This strongly suggests to try the following 
iteration scheme: 

h = n[F(f )], 

h =nnh)i 

(id 

/„ = K[F(f n -i)] 

Indeed, an iteration procedure of this type was successfully ap- 
plied in the so-called self-consistent field method ( Ostriker & Mark 

1 The fact of the surface of non-rotating stars is spherical allows us to spec- 
ify boundary conditions with one real number - the radius. Generally, in 
rotating objects we have a surface represented by some function of two 
variables which has to be determined. 

2 In differentially rotating stars the central density may be, but generally is 
not, the maximum density. 



ll968UHachisul fl986). We have introduced the canonical form to 
ensure that the first-order approximation is found in a correct order 
i.e. by using the first line of the sequence <1 It . When we go back 
to non-canonical form JS} the first line of <1 H takes the form: 

C-$ c -h(p 1 )=K( Po ) (12) 

From eq. 1121 above we can find the first-order deviation from 
sphericity. It seems impossible at first sight to avoid explicit inte- 
gration in eq. 1121 . In case of a general po, this is true. But let us 
look at equation J8j in case of vanishing centrifugal potential $ c , 
i.e. with no rotation: 

h(p)+K(p) = C (13) 

When we use a function which satisfies eq. 1 1 3i as zero-order ap- 
proximation: 

h(p )+1l(p ) = C (14) 
integration in eq. 1 1 2t can be easily eliminated: 

C-9 B -h{pi)=H{po) = Co-h{po) (15) 

Finally, our formula takes the form: 

h( Pl ) = h{p )-$ c + C-C (16) 

or simpler, using the enthalpy (h(po) = ho, h(pi) = hi): 

hi=h -$ c + C- C (17) 

Functions used as zero-order approximation (po or ho) are simply 
density and enthalpy distributions of non-rotating barotropic stars. 
In case of polytropic EOS, p — Kp 1 , these quantities are given by 
Lane-Emden functions. In more general case we have to find solu- 
tion of the ordinary differential equation of the hydrostatic equilib- 
rium. 

The only unanswered question is what are the value of con- 
stants C and Co 3 . It is an essential part of this work, so we decided 
to explain it in a separate section. 



4 ADJUSTING CONSTANTS 

When we try to find the enthalpy distribution using the formula 
J 171 we have to find the best zero-order function ho and the value 
of Co given by ho- An equivalent problem is to find the equation 
for which the function ho and the value Co are best zero-order ap- 
proximations - in this case we seek for C. We consider the latter 
case, as we have to find only one real number. Let us denote: 

AC = Co - C. (18) 

In our approximation, the equation for the first-order enthalpy dis- 
tribution hi, in terms of the initial spherical distribution ho and 
rotation law is, from 1171 : 

hi = ho - $ e - AC (19) 

where AC is still to be determined. 

The terms in eq. I19i behave as follows: 

(i) 'ho' is spherical enthalpy distribution, thus it is only a func- 
tion of the radius, has a maximum at the centre, and goes monoton- 
ically to zero, where usually is cut. However, from mathematical 
point of view, Lane-Emden functions extend beyond the first zero 
point with negative function values. 

3 For a given EOS the function po gives Co and vice versa. 
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Figure 1. Schematic graphs of the two first terms in eq. I19t . explaining the 
meaning of AC. The most general case is shown. In some particular cases 
C( ' andC' 2 ' may not exist- this depends on the both ho and & c . Vertical 
dot-dashed line shows that C' 1 ) = $ c (^?o). where Rq is radius of non- 
rotating star. The dashed curve fragment below the axis reflects ambiguity 
of the Lane-Emden function continuation to negative values, as described 
in the text. This figure prepared with use of eq. J40) for 2 = 0. 



(ii) '— $ c ' is a monotonically increasing function of distance 
from the rotation axis. It starts with zero at the rotation axis. It does 
not change the enthalpy along the axis of symmetry. The strongest 
enthalpy increase takes place along the equatorial plane. 

(iii) 'AC shifts the sum of positive functions ho and — <£> c 
down. 

At first sight, shifting down by AC seems not needed (i.e. one 
would adopt AC = ) , because we obtain correct qualitative be- 
haviour - the star is expanded along equator. But often it is enough 
to introduce slow rotation to get positive value of ho — <3> c (i.e. for 
our approximation to the enthalpy in this case) for any r > 0, z = 
i.e. equatorial radius becomes infinite. It leads directly to physically 
unacceptable results - infinite volume and mass. So the value AC 
plays a non-trivial role and has to be found. Fig. \\\ shows the be- 
haviour of all terms in eq. <19> along the equator of the star, where 
the rotation acts most strongly. Horizontal lines show points where 
enthalpy is cut for a given value of AC. 

We can distinguish some important values: 

(i) For AC < C (0) we obtain infinite radius of a star. These 
values obviously have to be rejected. 

(ii) For C (0) < AC < C (1) = $ c (Ro), where R is the radius 
of a zero-order density distribution, we get finite volume of a star, 
but we use extension of ho with negative values. This introduces 
some problems which we discuss later in the article, although the 
resulting enthalpy and density are positive and physically accept- 
able. 

(iii) For C (1) < AC < C (2) = h (r = 0, z = 0) we get a 
density distribution which is topologically equivalent to the ball. 

(iv) For C (2) < AC < C (3) = (ho - ^ c ) max we get toroidal 
density distribution. This case exists only if strong differential ro- 
tation is present. 

(v) For AC > C' 3 ' the star disappears. 

We expect to find the solution in the range C (0) < AC < C (2) 
because we are looking for finite- volume non- toroidal stars. 

One can try to find AC both analytically and numerically. To 
keep the algebraic form and the simplicity of the formula, we now 
concentrate on the former method. 



When we substitute the formula I19i into our basic equation 
^5) we get: 



h - $ c + C - C, + $ g (hi) + $c = C 



(20) 



In this formula we have made use of 1 1 81 . After obvious simplifi- 
cations, using 1 1 -ft and denoting IZ(po) = § g (ho) we have: 



* e (fci) = * 9 (fto) 



(21) 



This equality is true only if po = Pi- The same holds for the en- 
thalpy: 

hi = h . (22) 
Using formula 1 1 91 again we finally obtain: 

AC = -$ c (23) 

Left-hand side of eq. 1231 is constant, while the right-hand side is a 
function of distance from the rotation axis, monotonically decreas- 
ing from zero. This equality holds only in trivial case AC = and 
&c = with no rotation at all. In any other case 1231 cannot be 
fulfilled. So instead we try another possibility and require that 



AC: 



(24) 



where 'hat' denotes some mean value of the function $ c . We have 
chosen 

AC = -$; = - (Ji*Ro 3 } J $ c d 3 r (25) 

Vo 

Integration is taken over the entire volume Vo = |vri?o °f a non- 
rotating initial star with the radius Ro- This choice of AC gives 
good results. But using the mean value theorem: 



(26) 



where <& c is some value of $ c in the integration area, and taking 
in account monotonicity of the centrifugal potential < — <fr c < 
-$ c (Ro) we get: 

AC = -C < -$ c (Ro) (27) 

i.e. the value of AC is in the range AC < C (1) from Fig.Q] It 
forces us to use negative values of non-rotating enthalpy. Moreover, 
in case of polytr opic EOS with fractio nal polytropic index 4 Lane- 
Emden equation (KiDDenhahn & Wei g eri 19941) : 



z 2 dz 



< dw 

dz 



+ w = 



(28) 



has no real negative values, because of fractional power of negative 
term w n . But we can easily write equation, with solution identically 
equal to solution of Lane-Emden equation for w > 0, and real 
solution for to < e.g.: 



1 d 



i dw 







z 2 dz \ dz 
But, for example, solution of the following equation: 



1 d f 2 dw 
z 1 dz \ dz 



|n+l 



= 



(29) 



(30) 



is again identically equal to solution of Lane-Emden equation for 

4 Physically interesting cases like degenerate electron gas in non- 
relativistic case has fractional polytropic index n = 3/2. 
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w > 0, but differs from eq. 1291 for w < 0. Fortunately, difference 
between solution of eq. 1291 (Fig. Q below axis, dot-dashed) and 
eq. 1301 (Fig.Q solid) for w < is small if |iw| <C 1. Example 
from Fig. Q (for n = 1) is representative for other values of n. 
We will use form 1291 instead of the original Lane-Emden equation 
1281 for calculations in this article. 

To avoid problems with negative enthalpy we can put simply: 



AC = -$ c (7? ) 



(3D 



which is strictly boundary value C' 1 ' from Fig.Q The great ad- 
vantage of the eq. 13 1 1 is the possibility to analytically perform the 
integration of the centrifugal potential J2j for most often used forms 
of fi(r). In contrast, in formula 1251 . not only angular velocity pro- 
file 0, but also the centrifugal potential have to be analytically 
integrable function. In both cases 1251 13 1 i however, possibility of 
analytical integration depends on the form of fl(r). The value of 
AC from eq. 13 1 i also gives reasonable iso-density contours, cf. 
Fie. II ll andl8l but global accuracy is poor (Table0. 

As we noticed, the best value of AC in formula 1191 could be 
found numerically. For example, we can use virial theorem formula 



round numerically, for example, we 
for rotating stars (cf. lTassou l 2000) 



2E k -\E g \+3 / pcfr = 



(32) 



where E k and E g is the rotational kinetic energy and the gravita- 
tional energy, respectively. We define, so-called virial test parame- 
ter Z: 



2E k -\E 3 \+3(y-l)U 
\E a \ 



where we introduced internal energy 

1 



U = 



7- 



j3 

pa r 



(33) 



(34) 



Parameter Z is very common test of the global accuracy for rotating 
stars models. We may request that our enthalpy satisfy <32t . i.e. we 
choose AC from equation: 



Z (h - $ c - AC) = 



(35) 



We can find AC from equation eq. <35> numerically only. 

As it is shown on Fig. [2] we can find approximation of the 
rotating polytrope structure in form I19i satisfying virial theorem 
J32I up to accuracy limited only by numerical precision. 

Values of AC obtained with J 3 It . 1251 and from virial test 
<35> are compared on Fig. [3] Some of the global model properties 
are very sensitive to value of AC (cf. Figs|2|and|5j. 

Because virial test is unable to check accuracy of our model, 
we may also try to compare directly eq. i 191 with enthalpy distribu- 
tion from the numerical c alculations h nuln of e.g. lHachisul 1 1 98r3) . 
lEriguchi & Miille3ll^ 85) and find AC minimizing e.g. the follow- 
ing formula: 



[MAC) 



(36) 



This method however, requires numerical results (e.g. enthalpy dis- 
tribution) in machine-readable form. 



5 APPROXIMATE FORMULA ACCURACY 

In the above sections, we tried to be as general as possible. Now we 
give some examples, and test accuracy of approximation. 



AC 



Figure 2. Behaviour of virial test parameter Z vs AC for n = 3/2 poly- 
tropic model 141 1 with j -const angular velocity profile for = 1.5 and 
A = 0.2Rq. Proper choose of AC can give density distribution satisfying 
virial theorem to arbitrary accuracy. AC = 0.308 from eq. 13 1 1 is repre- 
sented by circle. Cross marks value AC = 0.269 given by formula 1251 . 
Virial theorem will be satisfied when we take AC, given by intersection of 
Z with horizontal axis i.e. AC = 0.189. 



In case of polytropic EOS p = Kp 1 the enthalpy is: 

Hp) 



x-y 7-1 

TP 



(37) 



7-1' 

Zero-order approximation of density (density of non-rotating poly- 
trope, |KiDDenhahn & Weigerl 1994) with n-th Lane-Emden func- 
tion w n is: 

4ttG 



po = Pc [w n (Ar)] n , j 
and our formula for density becomes: 



Pi = 



l/n 
pc Wn — 



nKj' 



($ c + AC) 



(38) 



(39) 



where AC is calculated from 1251 . 13 II or 1351 . 

In certain cases Lane-Emden functions are elementary func- 
tions as e.g. w%. In cases like this our formula may be expressed 
even by elementary functions. For example, for n = 1, K = 1/2, 
47rG=l, p c = l, ft(r) =fi /(l + r 2 /A 2 ) and AC from eq. J5TJ 
we get a simple formula: 



Pi(r,z) 



inVr 2 + z 2 1 n 2 A 2 r 2 1 n§A 2 7r 2 
y r 2 + z 2 + 2l + r 2 /A 2 2l + n 2 /A 2 



(40) 



Functions like this can easily be visualized on a 2D plot. Fig- 
ureQhas been made from the formula <40> while figures l8landfTT1 
from eq. I41i . 

Now we concentrate onn = 3/2 polytrope. In our calcula- 
tions and figures we will use 4-7rG = 1, p c = 1 and K = 2/5. 
Now formula I39i becomes: 

pi = 0„ - * c - AC) 3/2 (41) 
Iso-density contours of pi from J4H are presented on Fig. ll H and 
Fig.0 

To test accuracy of approximation we have calculated axis ra- 
tio, total energy, kinetic to gravitational energy ratio, and dimen- 
sionless angular momentum. Axis ratio is defined as usual as: 

R 

Axis Ratio = — — 

5 7 = 1 + ^ 



(42) 
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Figure 3. Dependence of AC on Qq, given by eq. 1311 (dotted), calculated 
from 1251 (solid) and given by virial theorem constrain 1351 (dashed). Both 
values estimated by ||5]^5j are below AC = -$ c (iio) i.e. C^ 1 ' from 
Fig.Q It shows, that continuation of Lane-Emden equation to negative val- 
ues is required for successful approximation of the rotating body structure. 
Density distribution was given by eq. 1411 with j'-const angular velocity 
with A = 0.2Rq. 



where R z is distance from centre to pole and R eq is equatorial 
radius. Total energy Etot- 



Etot = (E k + E 3 + U)/E 



is normalized by: 



E = (4ttG) 



2 M^_ 
J 2 



and dimensionless angular momentum is defined as: 



.2 
J = 



1 



J 



■P 1/3 

Hmax 



(43) 



(44) 



(45) 



4-kG M 10 /3 

where M and J are total mass and angular momentum, respec- 
tively; p m ax is maximum density. Quantities 1421 - 1451 are com- 
puted numerically from <4U . with given angular velocity fi(r) and 
chosen AC. 



Figure 4. E^ / \ E g | ratio as a function of the square of dimensionless an- 
gular momentum j 2 for our model 1411 with n = 3 /2, f2p = 1.5 and 
A = 0.2R(). Solid line represent numerical results of lEriguchi & ' Mullet 
Il985t) . We see that our formula behaviour is in good agreement with nu- 
merical results if E k /\E g \ <C 0.1. Results using AC from eq. 1251 are 
marked by crosses. Results satisfying virial theorem (AC from eq. 1351 ) 
are represented by diamonds. As it is apparent from figure above, AC has 
no influence on this relation, and can't improve accuracy of the formula 

en. 




5.1 Influence of AC 

We have made detailed comparison of our n = 3/2 model {41} 
with j-const rotation law and A — 0.2Ro (middle row of Fig. [8) 
for different values of AC with results of I Eriguchi & Mulle rfl985L 
Table lb). TableQshow our results for AC from eq. 1251 . Value 
of AC from eq. 13 II and corresponding virial test parameter Z is 
included here for comparison. Table [2] shows global properties of 
our approximation with AC equal to the solution of eq. 1351 . i.e. 
satisfying virial theorem. 

Direct comparis on of values from Table \\\ and table Ta- 
ble to Table lb of Erig uchi & Miillerl ll985t) may be difficult, 
becaus e our driving paramete r is central angular velocity Qq, 
while Eriguchi & Miiller 1 1985), following successful approach of 
lHachisul7l9 86), use axis ratio 1421 . More convenient in this case 
is comparison of figures p repared from data found in Table lb of 
lEriguchi & Mullerl jl985B and our tables. This is especially true, 
because axis ratio isn't well predicted by our formula (cf. Fig.|6|and 
Fig.0, while global properties (E t ot, j 2 , E k /\E g \, cf. Fig.|4||5} 
and virial test Z are in good agreement if Ek/\E g \ 0.1. 

Fig-0shows that our approximation is valid until Ek/\E g \ ~ 



Figure 5. Total energy versus j 2 . We can see great improvement of re- 
sults with AC from eq. marked by '+'. AC from eq. (©) 
gi ves in correct behavi our of the total energy. Solid line represent results 
of Eriguchi & Mullet] (1985J), cross and diamonds are result derived from 
our approximate formula 141 1 with AC from virial test 1351 and eq. 125 1 , 
respectively. Approximation satisfying virial equation gives results resem- 
bling numerical calculations. Parameters of model are given in caption of 
Fig.@ 

0.05, and begins to diverge from numerical results strongly for 
Ek/\E g \ > 0.1. Both values of AC t25!35> give similar behav- 
ior here. However, AC from virial test produces better results, and 
Ek/\E g \ values are more sensible for the strongest rotation. 

In contrast, total energy <43t is very sensitive to AC. Value of 
AC from eq. 1251 produces wrong result. E to t begin to increase for 
j 2 > 0.007, while numerical results give monotonically decreasing 
Etot- Use of AC from <35> instead, gives correct result, cf. Fig.|5| 

While global properties of our model are in good agreement 
with numerical results for Ek/\E g \ 0.1, axis ratio tends to be 
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n 


Axis Ratio 




j 2 


^ k 


E n 


+ E k + U 


Virial test Z 


Z 


AC 


AC 




\Eg\ 




E 


0.25 


1.01 


8.41 


x 10~ J 


0.004 


—2.50 


x 10 


O.Ol 


O.Ol 


O.Ol 


O.Ol 


0.50 


1.05 


3.38 


x 10 -4 


0.02 


— l.OO 


x 10~ 5 


0.03 


0.04 


0.03 


0.03 


0.75 


1.15 


7.64 


x 10~ 4 


0.04 


—2.20 


x 10~ 5 


0.07 


O.IO 


0.07 


0.08 


1.00 


1.19 


1.37 


x 10~ 3 


0.07 


—3.80 


x 10~ 5 


0.12 


0.17 


0.12 


0.14 


1.25 


1.31 


2.17 


x 10 -3 


O.IO 


—5.56 


x lO -5 


0.19 


0.28 


0.19 


0.21 


1.50 


1.46 


3.12 


x 10 -3 


0.15 


—7.38 


x 10~ 5 


0.27 


0.43 


0.27 


0.31 


1.75 


1.68 


4.40 


x 10 3 


0.20 


-8.93 


x 10~ 5 


0.36 


0.63 


0.37 


0.42 


2.00 


1.98 


5.90 


x 10~ 3 


0.26 


-9.89 


x 10~ 5 


0.47 


0.89 


0.48 


0.55 


2.25 


2.45 


7.55 


x 10~ 3 


0.32 


-l.Ol 


x lO- 4 


0.57 


1.22 


0.60 


0.69 


2.50 


3.30 


9.30 


x 10~ 3 


0.38 


-0.94 


x lO- 4 


0.67 


1.66 


0.75 


0.86 


2.75 


5.80 


1.10 


x 10~ 2 


0.44 


-0.84 


x lO- 4 


0.75 


2.22 


0.90 


1.04 



Table l. Properties of n = 3/2 polytropic model I41I with j -const rotation law and A = 0.2Rq. AC is calculated from eq. I25I and AC from eq. I ill 
Virial test in the latter case is labeled by Z. By little change from AC to AC one can notice significant improvement of virial test. In both cases virial test 
suggest strong deviation from equilibrium, especially for strong rotation. Actually, virial test is very sensitive to AC, cf. Fig.|2| We can require virial theorem 
to be satisfied, by use of another value for AC, solution of eq. 1351 . Results can be improved significantly this way - compare with Tablel2land Fig.l5l 




E k /IEgl 

Figure 6. Axis ratio versus / \ E g \ . We see that our formula (+ , o) under- 
estimates axis ratio. Choose of AC satisfying virial theorem (+) improves 
situation a bit. Solid line again is result of <Eriguchi & Miilleill 985). 




Figure 7. Axis ratio vs j 2 . Everything the same as in Fig. |S| AC from 
eq. 1331 (+) gives better approximation to axis ratio compared to formula 

E3 (o). 



underestimated, even for small values of j 2 . Fig.[6|and Fig.0show 
minor improvements when we use AC from virial test 1351 instead 
of mean value <25> . 

This subsection clearly show importance of constant value 
AC. Best results are produced with AC from eq. 1351 . therefore 
this value will be used in the next subsections to investigate influ- 
ence of differential rotation parameter A and type of rotation law 
on formula accuracy. 



5.2 Effects of differential rotation 

In addition to the results from previous subsection (j-const with 
A — 0.27?o) we have calculated properties of the almost rigidly 
(A = 2i? ) and extremely differentially (A = 0.02i? ) rotating 
model with the same rotation law. 

In all three cases we are able to find value of AC satisfy- 
ing eq. 1351 . However, this is not enough to find correct solution, 
because other parameters describing rotating body may be wrong. 
This is clearly shown on Fig. [5] where E^/Eg versus j 2 I45i is 
plotted for three cases of differential rotation. Apparent discrep- 



ancy for A = 2Ro exists. Both j-const and ?;-const angular veloc- 
ity profiles behaves as rigid rotation in this case. Thus we conclude 
that our formula is unable to predict correct structure in case of 
uniform rotation even if rotation is small. 

If rotation is concentrated near rotation axis, like in A = 
0.02iio case, our and numerical results are of the same order of 
magnitude. Quantitative agreement is achieved only for very small 
values of Qo- Let's note that in this case AC required by virial the- 
orem 1351 is slightly below zero (Table |4j- This example shows, 
that AC may also be negative. All three cases are summarized on 
Fig.|9l 

Results from this section show, that our formula is able to 
find correct structure of rotating body for differential rotation only. 
Range of application vary with differential rotation parameters, and 
best results are obtained in middle range i.e. A — 0.2Ro- With ex- 
tremal case (A = 0.02i?o) quality of our results is significantly 
degraded. 

In next subsection we examine, if this statement depends on 
rotation law. 
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Analytical approximation for rotating barotropes 1 



n 


Axis Ratio 


j 2 








_|_£; fc _|_[/ 
E 


Virial test \ Z\ 


AC 


0.25 


1.01 


8.30 x 10~ J 


0.004 


—2 


05 


x 10 


9 x 10 


0.004 


0.50 


1.05 


3.22 x 10 -4 


0.02 


-0 


98 


x 10 -5 


3 x 10 -4 


0.02 


0.75 


1.11 


6.86 x 10~ 4 


0.03 


— 2 


10 


x 10~ 5 


3 x 10~ 5 


0.04 


1.00 


1.20 


1.14 x 10~ 3 


0.06 


—3 


50 


x 10~ 5 


4 x 10 -5 


0.08 


1.25 


1.33 


1.65 x 10~ 3 


0.08 


—5 


14 


x 10~ 5 


6 x 10 -4 


0.13 


1.50 


1.51 


2.20 x 10~ 3 


0.11 


—6 


84 


x 10~ 5 


4 x 10 -4 


0.19 


1.75 


1.74 


2.77 x 10~ 3 


0.13 


-8 


55 


x 10~ 5 


3 x 10~ 5 


0.27 


2.00 


2.04 


3.32 x 10~ 3 


0.16 


-1 


00 


x 10- 4 


1 x 10" 3 


0.37 


2.25 


2.47 


3.89 x 10~ 3 


0.18 


-1 


18 


x 10~ 4 


9 x 10~ 4 


0.48 


2.50 


3.12 


4.40 x 10~ 3 


0.20 


-1 


33 


x 10- 4 


3 x lO" 4 


0.62 


2.75 


4.34 


4.93 x 10~ 2 


0.22 


-1 


50 


x lO" 4 


1 x 10" 6 


0.77 



Table 2. The same model as in Tabled] but no w AC is derived numerically from eq. 1351 . Virial test show accuracy of solution to eq. 1351 . Comparison 
with Tableland Table lb of lEriguch i & Muller <198.^ shows significant improvement of the total energy. Axis ratio is also closer to results of numerical 
calculations, and stability indicator £^/|.E g | isn't unreasonably high. See also Figs l4lT7l 




Figure 8. Examples of the density distributions given by our formula 1411 . Results for polytropic model with 7 = 5/3 and so-called j'-const lEriguchi'& Mullet 
1985) rotation law £2(r) = f2o/(l + r 2 /A 2 ) are presented. Upper row corresponds to diifferential rotation with A = 0.02Rq- The value of Slo increases 
from left to the right. Lower row shows behaviour of the almost rigidly rotating star with A = 2Rq. For the middle row A = 0.2i?o is adopted. Ro is the 
radius of non-rotating star. Values of the fl are, from left: upper row 75, 150, 200, 250; middle row 0.5, 1.0, 1.5, 2.0; bottom row 0.01, 0.02, 0.03, 0.035. 
Constant AC is calculated from eq. 1311 . 
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8 A. Odrzywotek 



Q 


Axis Ratio 


P 




B S 




K„ 4- E > A-U 




Virial test \Z\ 


AC 


O.Ol 


1.04 


6.10 x 10" 


-6 


1.9 x 10" 


4 


-1.83 x 10" 


7 


9 X 10~ 5 


0.01 


0.02 


1.20 


2.96 x 10" 


-5 


8.7 x 10" 


4 


-8.80 x 10" 


7 


1 x 10~ 4 


0.06 


0.03 


1.84 


1.26 x 10" 


-4 


2.9 x 10" 


-3 


-3.58 x 10" 


6 


5 x 10~ 4 


0.15 


0.035 


2.30 


4.49 x 10" 


-4 


7.0 x 10" 


-3 


-1.13 x 10~ 


5 


3 x 10~ 5 


0.24 



Table 3. Properties of our approximate sequence in case of j'-const rotation law with A = 2Ro, i.e. almost uniform rotation. In spite of fact that virial test is 
fulfilled with accuracy of order 10 -4 , comparison of data in this table with numerical results (cf. Fig.|9j clearly shows that our formula fails in case of rigid 
rotation. 



H 


Axis Ratio 


f 




Em 
W7\ 


E 


u 


Virial test \ Z\ 


AC 


25 


1.01 


1.45 x 10" 


-4 


0.02 


-0.44 x 10 


-5 


4 x 10" 4 


-0.01 


50 


1.05 


4.73 x 10" 


-4 


0.07 


-1.42 x 10 


-5 


3 x 10~ 4 


-0.02 


75 


1.12 


8.31 x 10" 


-4 


0.12 


-2.47 x 10 


-5 


4 x 10~ 4 


-0.03 


100 


1.23 


1.17 x 10" 


-3 


0.16 


-3.36 x 10 


-5 


2 x 10~ 4 


-0.02 


150 


1.55 


1.63 x 10" 


-3 


0.23 


-4.65 x 10 


-5 


4 x 10~ 4 


0.07 


200 


2.06 


1.96 x 10" 


-3 


0.28 


-5.50 x 10 


-5 


3 x 10~ 4 


0.24 


250 


2.95 


2.21 x 10" 


-3 


0.32 


-6.12 x 10 


-5 


1 x 10~ 4 


0.49 


300 


5.76 


2.42 x 10" 


-3 


0.34 


-6.64 x 10 


-5 


3 x 10~ 4 


0.83 



Table 4. Properties of sequence with j'-const rotation law for A = 0.02iJo 




y 

Figure 9. Stability indicator E k /E g vs j 2 for j'-const angular velocity 
(log — log plot) for three values of A. Quantitative agreement between our 
formula (symbols) and numerical results lEriauchi & Mullet 1985, lines) is 
achieved for A = 0.2iJo if Ek/Eg ^ 0.1. This case is presented as solid 
line and crosses. In case of A = 0.02iJo we have results of the same order, 
but they are identical only in case where rotation strength is very small. This 
case is presented by dashed line and diamonds. Formula fails (dotted line 
and circles) in case of rigid rotation. 



5.3 Rotation law effects 

In addition to previously described cases, we have calculated global 
properties of our model in case of n-const angular velocity profile, 
with parameter A = 0.2R (Tabled and A = 0.02i? (Table^. 
Results with A — 2Ro aren't presented, because they are similar 
to j-const case (cf. Table|3j, where both functions f2(r) behave as 
uniform rotation, and our formula fails in this case. 

Figures [Tol and fT2l show very good agreement of of the global 
physical quantities (Ek/E g , j 2 , Etot) with numerical results for 




Figure 10. E k / E g (j 2 ) for -u-const rotation law with A 



0.2R 



(dashed, cross) and A = 0.02i?o (solid, diamond), where lines refers to 
lEriguchi & Mulled Ii985t and symbols refers to our formula with AC from 
1351 . In this case we have good quantitative agreement with numerical re- 
sults in both cases up to E k /E g ~ 0.1. 



entire range of rotation strength covered by both methods. The most 
extreme case (A = 0.02) also behaves well. Axis ratio (Fig. 1131 
however, clearly distinguish between approximation and precise 
solution. Results are quantitatively correct only for small rotation 



parameters, e.g j 2 < 0.005 i.e. E k /E g < 0.1. 



6 DISCUSSION& CONCLUSIONS 

Comparison of the results obtained with our approximation for- 
mula (Fig.[TT]- Fig. [13] Table|2H5j with other lEriguchi & M uller 
Il985l Fig. 2-5, Fig. 9, Table 1 and 2) shows a correct qualitative 
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Analytical approximation for rotating barotropes 9 



n 


Axis Ratio 


P 




R. 
a k 

S.9 




Virial test \Z\ 


AC 


0.25 


1.05 


3.29 x 10" 


-4 


0.01 


-0.99 x 10~ 5 


9 X 10~ 5 


0.01 


0.50 


1.21 


1.41 x 10" 


-3 


0.05 


-4.26 x 10~ 5 


2 x 10~ 4 


0.05 


0.75 


1.63 


3.80 x 10" 


3 


0.11 


-1.12 x 10~ 4 


2 x 10~ 4 


0.13 


1.00 


2.31 


9.22 x 10" 


3 


0.21 


-2.52 x 10~ 4 


3 x 10~ 4 


0.27 


1.25 


5.26 


1.44 x 10" 


2 


0.31 


-3.87 x 10" 4 


5 x 10~ 4 


0.47 



Table 5. Properties of our model with f-const rotation law and A = 0.2Rq . 





Figure 11. Another example of density distributions given by our formula ED. The same as in Fig.|8]for so-called D-const rotation law i Eriguchi & Miillei 
1985) n(r) = n /(l + r/A). Values of U are, from left: in the upper (A = 0.02R ) row: 1.0, 3.0 , 5.0, 7.0; in the bottom (A = 0.2R ) row: 0.5, 0.75, 
1.0 and 1.25. 





Figure 12. Etot versus square of dimensionless angular momentum j 2 . 
Symbols description is the same as on previous figure, Fig. 1101 



Figure 13. Axis ratio versus j 2 for D-const rotation law. Symbols descrip- 
tion is the same as on previous figures, Fig. llOl and Fig. 1121 



behaviour for even the most simplified version of our approxima- 
tion formula for a wide range of parameters describing differential 
rotation and strength of rotation. This make our formula excellent 
tool for those who are interested in the structure of barotropic, 



differentially rotating stars, but do not need exact, high precision 
results. It can be applied for qualitative analysis of structure of 
rapidly rotating stellar cores (e.g. 'cusp' formation, degree of 
flattening, off-centre maximum density) with arbitrary rotation 
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10 A. Odrzywolek 



fl Axis Ratio j 2 4h B n+ E k+u virial test 121 AC 



1.0 1.02 9.81 x 10~ 5 0.004 -0.30 x 10~ 5 7 x 10~ 5 0.01 

2.0 1.09 3.94 x 10" 4 0.015 -1.20 x 10" 5 1 x 10" 4 0.04 

3.0 1.21 8.96 x 10~ 4 0.03 -2.75 x 10~ 5 7 x KT 4 0.09 

4.0 1.41 1.63 x 10~ 3 0.06 -5.00 x 10" 5 5 x 10~ 4 0.17 

5.0 1.73 2.65 x 10" 3 0.09 -8.09 x 10" 5 1 x 10" 4 0.27 

6.0 2.24 4.08 x 10" 3 0.13 -1.23 x 10" 4 3 x 10" 4 0.39 

7.0 3.03 6.05 x 10" 3 0.17 -1.77 x 10" 4 2 x 10" 4 0.55 

8.0 4.60 8.46 x 10~ 3 0.22 -2.38 x 10" 4 3 x 10~ 4 0.75 



Table 6. Properties of ii-const sequence for A = 0.02i?o. 



law, also for initial guess for numerical algorithms. It can also be 
used as an alternative for high-quality numerical results for use 
in a more convenient form as long as Ek/\E g \ <g 0.1 and we 
are interested mainly in global properties of differentially rotating 
objects. 
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